/** 
 * \file gdaxpy.cu
 * \author Kyle E. Niemeyer
 * \date 10/04/2011
 *
 * Based on "daxpy.f" from BLAS.
 *
 */

////////////////////////////////////////////////////////////////////////

/** gdcopy adds scalar multiple of one double precision vector to another. BLAS level one.
 * Uses unrolled loops for increments equal to one.
 *
 * \param[in]   n       vector dimension
 * \param[in]   dx       array, dimension (n-1)*|incx|+1
 * \param[in]   incx    increment between elements of x
 * \param[out]  dy       array, dimension (n-1)*|incy|+1
 * \param[in]   incy    increment between elements of y
 */
__device__ __inline__ void gdaxpy ( int n, double da, const double *dx, int incx, double *dy, int incy )
{
    
  if ( n <= 0 ) return;
  if ( da == 0.0 ) return;
  
  if ( ( incx == 1 ) && ( incy == 1 ) ) {
    
    // both increments equal to 1
    
    // clean-up loop
    
    // uint m = n % 4;
    // use faster bitwise operation instead
    uint m = n & 3;
    
    if ( m != 0 ) {
      for ( uint i = 0; i < m; ++i ) {
        dy[i] += da * dx[i];
      }
    }
    
    if ( n < 4 ) return;
    
    for ( uint i = m; i < n; i += 4 ) {
      dy[i] += da * dx[i];
      dy[i + 1] += da * dx[i + 1];
      dy[i + 2] += da * dx[i + 2];
      dy[i + 3] += da * dx[i + 3];
    }
    
  } else {
    
    // unequal increments or equal increments not equal to 1
    
    uint ix = 1;
    uint iy = 1;
    if ( incx < 0 ) ix = (-n + 1) * incx + 1;
    if ( incy < 0 ) iy = (-n + 1) * incy + 1;
    
    for ( uint i = 0; i < n; ++i ) {
      dy[iy] += da * dx[ix];
      ix += incx;
      iy += incy;
    }
    
  } // end increment if
  
} // end gdaxpy